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Abstract 

When the frequencies of the elastic and pendular oscillations of an elastic pendu- 
lum or swinging spring are in the ratio two-to-one, there is a regular exchange of 
energy between the two modes of oscillation. We refer to this phenomenon as pulsa- 
tion. Between the horizontal excursions, or pulses, the spring undergoes a change of 
azimuth which we call the precession angle. The pulsation and stepwise precession 
are the characteristic features of the dynamics of the swinging spring. 

The modulation equations for the small-amplitude resonant motion of the system 
are the well-known three- wave equations. We use Hamiltonian reduction to deter- 
mine a complete analytical solution. The amplitudes and phases are expressed in 
terms of both Weierstrass and Jacobi elliptic functions. The strength of the pul- 
sation may be computed from the invariants of the equations. Several analytical 
formulas are found for the precession angle. 

We deduce simplified approximate expressions, in terms of elementary functions, 
for the pulsation amplitude and precession angle and demonstrate their high accu- 
racy by numerical experiments. Thus, for given initial conditions, we can describe 
the envelope dynamics without solving the equations. Conversely, given the param- 
eters which determine the envelope, we can specify initial conditions which, to a 
high level of accuracy, yield this envelope. 
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1 Introduction 



The present work is concerned witli ttie tliree-dimensional motion of the elas- 
tic pendulum or swinging spring in the case of resonance. It continues the 
investigation described in previous studies by Lynch [13] and by Holm and 
Lynch [9]. In particular, the exchange of energy between quasi- vertical and 
quasi-horizontal oscillations and the stepwise precession of the swing plane 
are investigated. 

When the ratio of the normal mode frequencies of the spring is 2:1, a reso- 
nance occurs, in which energy is transferred periodically between vertical and 
horizontal oscillations. The first study of this resonance was that of Vitt and 
Gorelik [17]. We refer to the regular exchange phenomenon as pulsation. The 
motion has two distinct characteristic times, that of the fast oscillations and 
that of the slow pulsation envelope. As the oscillations change from horizon- 
tal to vertical and back again, it is observed that each horizontal excursion 
or pulse is in a different direction. We call this change in azimuth the preces- 
sion angle. The motion thus has three components: oscillation (fast), pulsation 
(slow) and precession (slow), closely analogous to the rotation (fast), nutation 
(slow) and precession (slow) of a spinning top [3] . ^ 

We consider two complementary questions, one direct and one inverse: 

Question 1. Given initial conditions, can we describe the envelope dy- 
namics without solving the equations? 

Question 2. Given the parameters which determine the envelope, can we 
specify initial conditions which yield this envelope? 

We provide a complete answer to Question 1. Analytical expressions are de- 
rived for the pulsation amplitude, precession angle and period in terms of the 
invariants of the motion. We also develop accurate approximate expressions 
for the pulsation amplitude and precession angle. Thus, the envelope dynam- 
ics may be deduced from the initial conditions. Question 2 is more recondite, 
but we can give a positive answer for the physically interesting case of strong 
pulsation. We derive approximate expressions for the angular momentum and 
Hamiltonian in terms of the pulsation amplitude and precession angle. Initial 
conditions can then be determined which yield the desired envelope to a good 
level of approximation. 

We briefly outline the contents of the paper below. When the amplitude is 
small, the Lagrangian may be approximated to cubic order. When it is aver- 

http://www.maths.tcd.ie/~houghton (Conor Houghton). 

^ A Java Applet illustrating the pulsation of the swinging spring may be found at 
|http : / / www ■ maths . t cd . ie/ ~plynch/SwingingSpring/ SS-Home-Page . html 
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aged over the fast oscillation time, a set of equations for the envelope ampli- 
tudes is obtained. These modulation equations, the three-wave equations, are 
presented in §2. They are found to have three independent constants of mo- 
tion and are therefore completely integrable. Small-amplitude perturbations 
about steady-state solutions are studied in §3, and a crude estimate of the 
precession angle is obtained. The general solution of the three-wave equations 
for finite-amplitude motions is derived in §4. The amplitudes are expressed in 
terms of elliptic functions and the phase angles as elliptic integrals. Analytical 
expressions for the stepwise precession of the swing-plane are then derived. 



Recently, Dullin et al. [5] constructed a canonical transformation in which 
the angle of the swing plane is a coordinate in an action-angle system. They 
showed that the precession angle is one of the two rotation numbers of the 
invariant tori of the integrable system. They obtained a simple equation for 
the precession angle by approximating an elliptic integral. They proved an- 
alytically that the resonant swinging spring has monodromy and concluded 
that the system provides a clear physical demonstration of this phenomenon. 



Several approximate expressions for the precession angle, involving only ele- 
mentary functions, are obtained in §5. One of these is equivalent to the for- 
mula reported in [5]. The approximate solutions are compared to the values 
obtained from the analytical expression, and are found to give remarkably ac- 
curate results. The intensity of the pulsation envelope is determined by solving 
a cubic equation whose coefficients are defined by the invariants. Thus, the 
direct question is fully answered, in the affirmative. 



To answer the inverse question, we assume the pulsation amplitude and pre- 
cession angle are given and derive expressions for the invariants. From these, 
appropriate initial conditions are easily determined. The expression for the 
precession angle is easily inverted. To obtain an invertible expression for the 
pulsation amplitude, we approximate the cubic by a quadratic, and obtain in 
§6 simple approximate expressions for the angular momentum and Hamilto- 
nian. These approximations may be used to control the envelope dynamics by 
an appropriate choice of initial conditions. 



In the concluding section, §7, we present a schematic diagram which shows the 
qualitative dependence of the envelope motion on the values of the invariants. 
This allows us to determine, at a glance, the general character of the solution 
for given values of the constants of motion. Several important special solutions 
are indicated on the diagram. 
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Fig. 1. The swinging spring. Cartesian coordinates are used, with the origin at the 
point of stable equihbrium of the bob. The pivot is at point (0, 0,£). 

2 The Dynamical Equations 

The physical system under investigation is an elastic pendulum, or swinging 
spring, consisting of a heavy mass suspended from a fixed point by a light 
spring and moving under gravity, g (Fig. 1). We assume an unstretched length 
£o, length I at equilibrium, spring constant k and mass m. The Lagrangian, 
approximated to cubic order in the amplitudes, is 

L = 1 (±2 + ^2 ^ i2) _ 1 [a;|(a;2 + y"^) + ^2 ^2] ^ 1^(^2 ^ ^2)^ ^ 

where x, y and z are Cartesian coordinates centered at the point of equilibrium, 
— \J qI^ is the frequency of linear pendular motion, ujz = ^k/m is the 
frequency of its elastic oscillations and A = £o^|/^^- The equations of motion 
in cartesian, spherical and cylindrical coordinates may be found in [13]. There 
are two constants of the motion, the total energy and the angular momentum 
about the vertical, and the system is not integrable. Its chaotic motions have 
been studied by many authors (see Refs. in [14]). 

2.1 The Time-averaged Equations 

We confine attention to the resonant case uiz = '^^r and apply the averaged 
Lagrangian technique. The solution is assumed to be of the form 
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x^^{a{t)exp{iuRt)} , (2) 
y^^{b{t)cxp{tu;Rt)}, (3) 
z = K{c(t) exp{2iujRt)} . (4) 

The coefficients a{t), b{t) and c{t) are assumed to vary on a time scale wfiicfi is 
mucfi longer than the time-scale of the oscillations, r = I/cor. The Lagrangian 
is averaged over this time, yielding 

(L) = lujRh{da* + bb* + 2cc*} + ^{K{a'^ + b'^)c*} 



where k, = A/(4cjr). The resulting Euler- Lagrange equations are the modula- 
tion equations for the envelope dynamics: 



ia = K,a*c, (5) 
ib — Kb*c, (6) 
ic^lK{a^ + b^). (7) 

2.2 The three-wave equations 

We now transform to new variables 

A^lK{a + ib), B^\K{a-ib), C^kc. (8) 

Then the equations for the envelope dynamics take the form 

iA = B*C, (9) 
iB = CA*, (10) 
iC = AB, (11) 

These three equations for the slowly- varying complex amplitudes A, B and 
C are the three-wave equations. The relevance of these equations in various 
physical contexts is discussed in [9]. They govern quadratic wave resonance 
in fluids and plasmas. Their application to resonant Rossby wave triads is 
considered in [15]. In Appendix A, we show that they are a special case of 
the Nahm equations which are used to construct soliton solutions in certain 
particle field theories. For further references to the three-wave equations and 
a discussion of their properties see [2] . 

The three- wave equations conserve the following three quantities: 
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H^^{ABC* + A*B*C)=U{ABC*}, (12) 
A^=|A|2 + |5|2^-2|C|^ (13) 
J^\A\^-\Bf. (14) 

The equations are completely integrable. The following positive-definite com- 
binations of and J are physically significant: 

N+ = ^{N + J)^ + |C|2 , N_ = l{N-J)^ \B\^ + |Cp . 



These combinations are known as the Manley-Rowe relations. Together with 
the Hamiltonian H, they provide three independent constants of the motion. 
We note that H is invariant under the symmetry transformations 



(A,S,C)^(Ae*^,5e-*^,C), 
{A, B, C) ^ {Ae''', B, Ce^>^) , 
{A,B,C)^{A,Be''',Ce'''). 



(15) 
(16) 
(17) 



These symmetries are associated, via Noether's theorem, with the three in- 
variants {J, A^+, N^}. Any two of the transformations generate the third. This 
reflects the inter-dependence of J, and N_. 



2.3 Reduction of the system 



To reduce the system for if 0, we express the amplitudes in polar form: 



A=|^|exp(iO, (18) 
5 = |B|exp(ir/) , (19) 
C = |C|exp(i7). (20) 

In general, the phases of A, B and C are not periodic. However, Q = J — iC + v) 
is periodic. The Hamiltonian may be written 



H ^\A\\B\\C\ cos C- 



The amplitude \C\ will be obtained in closed form in terms of elliptic functions. 
Once \C\ is known, \A\ and \B\ follow immediately from the Manley-Rowe 
relations 

\A\^^N+-\C\\ \B\^^N^-\C\K 
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N+ C\ 

Fig. 2. Polynomial $o as a function of |Cp. 

The phases ^ and r] may now be determined. Using the three- wave equations 
(9)-(ll) together with equations (18)-(20), we find 



H 



UP ' 



77 = 



ISP ' 



(21) 



so that ^ and 77 can be obtained by quadratures. Finally, C, is determined 
unambiguously by 



d\ci 

dt 



= -2ii"tanC and = |^| |5| |C| cosC ■ 



(22) 



It also follows from (21) and (22) that 



lAP \B\ 



ICP 



(23) 



The phase of C follows immediately, 7 = C+^+C) ^^'^ we can then reconstruct 
the complete solution using (18)-(20). 



2.4 The equation for \C\ 



Prom the equation (11) for (7, and its complex conjugate we get 

d\C\^ 



dt 



= 2Q{ABC*} . 



(24) 



Using the definition of the Hamiltonian, it follows that 

\A\^\B\^\Cf ^H^+ ['^{ABC*}]^ . 
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Applying this to the square of (24) and using the definitions of the Manley- 
Rowe quantities immediately yields an equation for |Cp alone: 

/ci|C|2y 



dt 



{N+ - \C\'){N^ - \C\')\C\' - H' . (25) 



We define the cubic polynomial ^odC'P) (plotted in Fig. 2) by 

$o(|C|2) = {N+ - |C|')(Ar_ - |C|')|C|2 . (26) 

Then the right hand side of (25) may be written 4[$o(|C'P) — H'^]- For small 
H^, this cubic has three positive real roots. If these roots, in descending order 
of magnitude, are denoted C^, C| and C|, it follows that 

0<C^<Ci<N_<lN <N+<Cf<N. (27) 

(We have assumed without loss of generality that J > 0). In the case of equal- 
ity of roots, the solution may be obtained in terms of elementary functions. 
We assume in general that this is not so and solve for |Cp in terms of eUip- 
tic functions. However, before doing this, we investigate perturbation motion 
about steady solutions. 



3 Small- Amplitude Modulation of Steady States 



We consider the case where the variations of the amplitudes about their mean 
values are small. This enables us to make additional approximations and derive 
simple estimates of the pulsation period and rate of precession. From these 
two quantities, the precession angle follows immediately. 



3.1 Steady State Motion 

We first consider solutions for which the amplitudes \A\, \B\ and \C\ are con- 
stant. The simplest cases are where the phases are also constant; then the 
three-wave equations become 

B*C = CA* = AS = , 
which give three particular solutions 
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(i) 


A 


= Aq, 


B = 


C = 





(n) 


B 




C = 


A = 





in) 


C 


— Co , 


A^ 


B = 






The first two solutions correspond to conical motions: the bob moves in a circle, 
clockwise or anti-clockwise, while the spring traces out a cone. These solutions 
are stable to small perturbations. The third particular case represents purely 
vertical oscillations; this motion is unstable [13]. 

More generally, from (22), constancy of the amplitudes imphes ( — ^—{C+v) ~ 
so that H — \A\\B\\C\, and the three- wave equations become 



-m=\B\\c\ 

-\B\f,= \C\\A\ (28) 
-\C\i=\A\\B\. 

Differentiating (25), a simple algebraic manipulation yields 

|C|2 = Co^ = ^ (2N - VN^ + 3 J2) (29) 



The other amplitudes are given by 



\Af 




1 r 

6 


\B\' 


= Bl^ 


1 r 

6 



{N + 3 J) + ^^WT^ 



{N - 3 J) + VN^ + 3 J2 



These solutions are the elliptic-parabolic modes (EP-modes) studied by Lynch 
[13]. The precession rate is given hy ^l = 4> — |(C [9]- From (28) it follows 
that 



For J = we have planar motion with 



These are the cup-hke and cap-hke solutions of Vitt and Gorehk [17]. 
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3.2 Perturbation about Elliptic- Parabolic Motion 



We consider small deviations about the steady EP-modc solutions. We write 
Id 2 = Ci + e where is given by (29) and |e| < Cl Then, if (25) is 
differentiated and nonlinear terms in e are omitted, we obtain 

fl + (2Vm + 3J^)e^0. (31) 



The solution is e{t) — e(0) cosa;pf, an oscillation about Cq with the pulsation 
frequency 

OOP = V2</m + 3 J2 . (32) 



For the EP-modes, the horizontal projection is an ellipse precessing at a con- 
stant rate ft. The perturbation is a pulsating motion, with sinusoidal time 
variation, in which the major and minor axes of the elhpse alternately expand 
and contract with period Tp = 27i/ujp. The area of the ellipse is proportional 
to J and remains constant [9]. It is straightforward to derive expressions in 
terms of elementary functions for the remaining amplitudes and the phases, 
but they are not required to determine the precession angle. 

We note from (32) that y/2N < ujp < 2\/2N. From the precession rate and 
the pulsation frequency, the precession angle follows immediately: 

A0 = nTp . 



Using (29), (30) and (32), this gives us 



2H Up 3 \V8H, 



2N - VA^^ + 3J2" 

+ 3 J2 



(33) 



For small angular momentum J <^ N, the term in square brackets is close to 
^ and 
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4 Analytical Solution of the Three-wave Equations 

4.1 Solution in Weierstrass Elliptic Functions 

We now derive an explicit analytical solution for |Cp, valid for finite ampli- 
tudes. The solutions for |Ap and follow immediately from the Manley- 
Rowe relations. Then (21) are integrated for the phases. The integrals turn out 
to be similar to those occurring for the spherical pendulum, so the approach 
of Whittaker [18] applies. The required properties of the Weierstrass elliptic 
functions are given in Whittaker and Watson [19], Ch. 20 and in Lawden [12], 
Ch. 6 (see also Abramowitz and Stegun [1], Gradshteyn and Ryzhik [7] and 
Byrd and Friedman [4]). 

4.1.1 Solution for the amplitudes 

The quadratic term on the right of (25) is removed by a simple transformation 
u = [C[7-/V - 1/3 and r = ^/Nt. Then we obtains 



This is the standard form of the equation for Weierstrass elliptic functions. 
The constants 512 and g^, called the invariants, are given by 



For small H^, the discriminant A = — 27 is positive and the three roots 
are real. This is the case of physical interest, and we assume the roots of the 
cubic are ordered so that ei > 62 > 63. Note that 61 + 62 + 63 = 0. The general 
solution of (35) is 

u = p{t + a) 

where a is an arbitrary (complex) constant. The function p{z) is defined by 




4{u - ei){u - 62) (m - 63) . 



(35) 





(36) 
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3{z} -2 



3{z} -2 



Fig. 3. Weierstrass's p-function with half-periods toi = 2, = i on the domain 

{z = x + iy : X £ [—2, +6], y G [—1, +1]. Left panel: real part; right panel: imaginary 
part. Values foi z = x — los are plotted as heavy lines. Calculations are based on 
(36). The function has double poles at 2; = 2ma;i + 2nw2- 

where the summation is over all integral m, n except m = n = 0. It has poles 
on the real line and is doubly periodic: p{z + 2moui + 2nuj2) — p{z) for all 
integers m and n. The difficult problem of determining uji and UJ2 from the 
invariants is discussed in §21.73 of [19]. The quantity u)^ is defined by requiring 
(jji + UJ2 + = 0. It may be shown that 

p{uji) = ei p{uj2) = 62 pM = 63 . (37) 



In the present case, uoi is real and CJ3 is pure imaginary (explicit expressions are 
given below). On the real line, p{z) is real, with values in the range [ci, +00). 
On the line z — uJz-\- x \i takes real values in the interval [ea, 62]. Moreover, as 
z varies along the edge of the rectangle from to uj-i to uji +a;3(= —UJ2) to uj^ 
to 0, p{z) is real and decreases monotonically from +00 to ci to 62 to 63 to 
—00. To satisfy the initial conditions, we choose a = uj^ — Tq, where Tq is real 
and may be taken as zero by a suitable choice of time origin. Then p(r + ui^) 
is real and oscillates between 63 and €2- The solution for the amplitude is 

|C|2 = 7V[| + p(T + a;3)l ■ (38) 



The behaviour of the Weierstrass p-function is shown in Figure 3. For general 
z it takes complex values. On the line ^{z} = ur^ the function is real with 
periodic oscillations, as indicated by the heavy lines at the front of the figure. 
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4- 1-2 Solution for the phase angles 
Weierstrass's zeta function is defined by 

f^^-piz), Jim [C(^) - z-^] = . (39) 

It is quasi-periodic in the sense that 

C(^ + 2a;i)=C(^) + 2C(u;i). (40) 

We note that ({z) is an odd function of z and will use the relation 

- (^2C{'^i) = l^ri . (41) 

The sigma function is defined by 

^loga{z) = az), lim^ = l. (42) 
az ^—'■0 z 

It is also quasi-periodic, such that 

a(z + 2u;i) = - ex.p[2((u;i)(z + ui)] a(z) . (43) 

Three other sigma functions may be defined. The relationship between the 
sigma functions and the Weierstrass p-function is similar to that between the 
theta functions and the Jacobi elliptic functions. We will not require these. 
We will require the identity 

^'^"^ -az-a)-az + a) + 2aa) (44) 



p{z) - p{a) 



(this follows from a consideration of the poles and zeros of the functions on 
each side). 



The solution (38) leads to a solution for \A 



2. 



2 N + 3J 



\^\ -—^ Npir + us). (45) 



Substituting in the first of (21) we have 



dr 6Np{T + 003) -{N + 3 J) 
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Now we introduce auxiliary constants k± defined by 



, , iV + 3J , , iV-3J 



6N 



6N 



Using (35), it follows that 



4IP 
7V3 



We must determine which sign for the derivatives should be chosen. Prom (27) 
the following sequence of inequalities holds: 



— i < 63 < 62 < e_ < i < e+ < ei < 



(46) 



Since 62 < e_ < e+ < Ci, it follows that k± lie on the line between Ui and 
LUi+uj-i, which determines the sign of the derivatives to be p'(k±) = 2iH/N^/^, 
a positive imaginary number. The equation for ^ thus becomes 



dr \2iJ p{r + U3) - p{k+) 



Using (44) this may be expressed in terms of zeta functions and using (42) it 
may be integrated immediately to yield 



2i 



log 



a{T + a;3 + «;+) 



+ 2C(«^+)T 



(47) 



A similar expression holds for r] — rjo with replacing Thus we obtain 
the expression for the azimuthal angle 0: 



cf>- 



cr(T + 003 — K+) (7(r + a;3 + k_) 
a{T + u;3 + K+) (t(t + a;3 - K_) _ 



This is the solution for the azimuth as a function of time. Using the quasi- 
periodic properties (40) and (43), the change in when r varies by 2u!i may 
be computed: 



= -iuJi{C{K+) - C(«-)) + <(<^i)(k+ - K_) 



(48) 
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This is the desired analytical expression for the pulsation angle. ^ 

We note two obvious special cases of (48). When J = we have 

yielding a zero result for A0. When H — 0, we have = cui and = —002, 

so 

A0 = -^[a;iC(a;2) - a;2C(a;i)] = ^ , (49) 



where we have used (41). These two special cases intersect in the homoclinic 
orbit (with J = H = 0) which has an infinite transition time. 



4-2 Solution in Jacobi Elliptic Functions 



While (48) is the analytical solution, it is not immediately obvious how nu- 
merical information may be extracted from it. The quantities on the right side 
are all computable in principle, but at the expense of considerable effort. It is 
therefore useful to seek an alternative expression, in terms of Jacobi elliptic 
functions. 



4.2.1 Solution for the amplitudes 

Recall that with the transformation u — \C\'^/N —1/3 and r = y/Nt, (25) was 
transformed to (35), which we write again for convenience: 

(^] ^A(u-e,){u-e2){u-es). (50) 



For solutions of physical interest, is sufficiently small that the three roots 
of the cubic are real. Defining the quantities 

^2 _ / _2 3 \ ^2 _ _ ^\ ^ 

Vei - 63/ 



a further transformation. 



w — 



u- Cs 



62 - 63 



^ The apparent discrepaney with the result of Whittakcr for the spherical pendulum 
(p. 106 in [18]) arises from our choice of convention that '^{uj^/loi} > 0. Our result 
is consistent with the rigid body formula (7.3.24) in Lawden [12], who adopts the 
same convention as we do. 
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brings equation (50) to the standard form 



The solution is w = sn(s — sq), or 
= 63 + (62 - 63) sn^(s - So) 



where sq is arbitrary. The Jacobi eUiptic function sn s has period AK, where 

/dw 
I , (52) 



so sn^(s — So) has period 2K. For definiteness, we set so = 0, which means 
choosing the origin of time where the solution has a minimum: 

= Cl + {Cl - Cl) sn'(z/ViVt) . (53) 



Clearly, \C\ oscillates between C3 and C2 with physical period 

T^2K/vVn. (54) 

The remaining amphtudes, \A\ and \B\, follow from the Manley-Rowe rela- 
tions: 

\A\^ = N+-\C\^ , \B\^ = N_-\C\^ . 



They have the same period as \C\ but vary in anti-phase with it and in phase 
with each other. We denote the minimum and maximum values of |yl| by ^3 
and ^42, and similarly for \B\. Thus 

N^ = Al + Cl = Al + Cl, N_ = Bl + Ci = Bl + Cl 



The initial values of the amphtudes (for so = 0) are 

1^(0)1=^2, 15(0)1 =i?2, \C{0)\ = C,. 



We note here an important scaling invariance of the three-wave equations. If 
the amplitudes are magnified by a constant factor and the time is contracted 
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by the same factor, the form of the equations (9)-(ll) is unchanged. Thus, the 
period of the modulation envelope motion varies inversely with its amplitude. 
The overall scale may be measured by \fN and the inverse dependence of T 
on this is seen in (54). 

The solutions (38) and (53) must be equivalent. This follows from identities 
relating Weierstrass and Jacobi elliptic functions. The complimentary modulus 
is defined as k' = Vl — /c^, and we write K' = K{k'). The parameters are 
related by 



, /e2 - 63 , lei - 62 K iK' 



) — W 1 — I ) ^6 — I 

61-63 \ ex- 63 V61 - 63 V61 - 63 



([7], p. 919). Then we have 

ei - 63 



p{z) = 63 + 



sn2(^ei - e^z) 



But the Jacobi function sn(s + iK'^ is given in terms of its value on the real 
line by 

1 

sn(s + iK') 



ksns 

and the equivalence between the two forms of solution follows immediately. 
4.2.2 Solution for the phase angles 

It remains to determine the phases. Integration of (21) furnishes the angles ^ 
and rj. We define 



2 Cj - Cj 62 - 63 . H H/N^l^ 



N+ - C| 6+ - 63 ' ^ v^/NAl V61 - 63(6+ - 63) ' 

2 CI - CI 62-63 . H H/N^l^ 

7- = T7 7^ = Z ^ 



N_ - CI 6_ - 63 ' v^fNBl V61 - 63(6- - 63) ■ 

It follows from (46) that A;^ < 7^ < 7^ < 1. We may now write (21) in the 
form 

ds 1 — 7+ sn^s ' ds 1 — 7?. sn^s ' 
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The right sides are the integrands occuring in Legendre's elhptic integral of 
the third kind ([1], p. 590). They may be put in standard algebraic form by 
defining x = sns. Writing 



dx 



sn^s ^ (l-aa;2)J(l-a;2)(l-fc2a;2) 



the solution for ^ becomes 
e-eo = -A+n(s,7^,A;) 



(56) 



There is an analogous solution for rj. The changes in ^ and r] over a half period 
s e [0, K] are 

|Ae = -A+n(72 , k) , lArj = -A_n(72 , k) , 



where the complete elliptic integral is defined as Il{a,k) = Il{K,a,k). The 
azimuthal angle of the pendulum is — ^{C — v)- Thus, the change in the 
azimuth over a full pulsation period is 



A0 = -(A+n(72,fc)-A_n(72,fc 



(57) 



In Appendix B, an alternative formula (B.4) is derived from the expression 
(57), which is structurally similar to (48) obtained above. Using this formula, 
the limiting cases A(j) = for J ~ and A0 = 7r/2 for H — are again 
derived, in agreement with (49). 



5 Approximate Formulas for the Precession Angle 



We have derived an exact analytical expression for the precession angle, involv- 
ing elliptic integrals. It is of interest to obtain more convenient approximate 
formulas, involving only elementary functions. It might be expected that the 
easiest way to do this would be to approximate (57) directly. However, it turns 
out that it is easier, and more transparent, to return to the differential equa- 
tions governing the the system, use them to write down an integral for the 
precession angle and approximate this integral. 

The precession angle is (f) — ^{^ — rj). Combining the two components of (21), 
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we obtain 

dt ~ 2\Am\' ■ ^^^^ 



Eq. (24) may be written 

= ±2^|A|2|S|2|C|2-i/2 _ (59) 

CLv 

Taking the quotient of these two equations, we get 

dd) JH 



(60) 



The pulsation of the amphtude \C\ occurs between C2 and C3, where C| and 
C| are the two smallest zeros of the polynomial 

$ = |>insnc|2 -H^^ |C|'(|C|' - A^+)(|C|' - N_) - . (61) 



It is also useful to write ^ — ^0 — H"^ where 

$0 = \A\^\B\^\C\^ = |C|2(|C|2 - N+){\C\^ - N_) 

is as defined by (26) and illustrated in Fig. 2. The two signs in the differential 
equation (60) correspond to phase changes during alternate half-cycles of the 
pulsation. The integral of (60) over a full cycle may be written formally: 

2 J \A\^\B\^^/^ ^ ' 

c| I II I ^ 

It is convenient to change the integration limits; to do this, we consider (62) 
as an integral over the complex Z- plane, where Z = |C|2 on the positive real 
axis. This gives 

, , JH r dZ 

^<P^-r 63 

4 J {Z -N+){Z -N_)J^Z) 



The contour ^1 encircles C| and C| and the square root in the integrand has 
two branch cuts, one from C| to C| and the other from Cl to +00. This is 
illustrated in Fig. 4. In addition to the three branch points, the integrand has 
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■^1 



• •• 



N+ 



C\ 



Fig. 4. Contours and '^2 in the Z-plane. 

two simple poles at Z = and Z = N_. In fact, the residues at these two 
poles sum up to zero: 



Res(iV+) = -Res(iV_) = 



-iJ 



1 



4 N+-N^ 



"4 



(64) 



Furthermore, the integrand goes to zero sufficiently fast as |Z| — >■ 00 that we 
can replace the contour ^1 by ^2 (Fig. 4). Returning to the original integral 
(62), this corresponds to a change of integration range to 

A(p= / ' ' ^ . (65) 



This interval is more convenient than the previous one because the integrand 
is small everywhere except near the lower limit of integration and because 
the point of inflection in $ can cause difficulties when approximating $ near 
C|. Since the integrand is dominated by its behaviour near the obvious 
approach would be to find a quadratic which approximates $ near this point. 
As Cl is the root of a cubic, it can be written in terms of iJ, J and N ^ but 
this expression is cumbersome and does not yield a convenient approximation. 
It is simpler to consider the behaviour of $ at — {N + J)/2. This point is 
close to Cf because H'^ must be small compared to N for the periodic motion 
to exist. ^ 

Having decided to approximate at rather than Cf, the next step is to 
approximate $0 — Z{Z — N_^)(Z — N_) by a quadratic with a root at N^: 

^>^ = Z^{Z-N+){Z-Z^). (66) 



It is possible to perform the resulting approximate integral. However, the 
solution is complicated unless Zi = N_ (see Appendix C). Thus, we consider 

'^o^Zo{Z-N+){Z-N_). (67) 

^ It can be shown easily that the maximum allowed value of H'^/N^ is Hqq = 
1/54 0.0185 and occurs for J = 0. 
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The quadratic and cubic $o both vanish at Z = and Z — A^_. They 
are also equal when Z — Zq. We consider two choices of Zq. 



First, we choose Zq to be the mean oi Z — and Z — that is Zq — N/2. 
The integral (65) becomes 

^ , JH 7 dZ 

A0 = — / , . (68 

2 J (7. - N.\(7, - N \ J-<DJ7,\ - m 



[Z-N+){Z-N_)^^o{Z)-m 

where Z^ is the larger root of '^q — H'^ — 0. Defining a — 2Z — N , we get 

A , T 2V2JHda , , 

A(^= / / =, 69 

i_ (a2- j2)^a2-(j2 + 8i/ViV) 

where (7+ = y^J^ + SH'^/N. This may be integrated analytically ([6], p. 72) to 
give 



2-(j2+8H2/-^) 



+ CXD 



tan 



\Vnj 



(70) 



Noting that tan ^ x = | — tan -"^ |, the phase change over a full cycle is 



This elegant approximate formula for the pulsation angle was reported by 
DuUin, Giacobbe and Cushman [5] , although they did not give the factor -s/S 
explicitly, and we refer to it as the DGC formula. Numerical experiments 
indicate that it is of high accuracy throughout the accessible domain. 

An alternative choice of quadratic approximation requires 'I'o and $o to have 
equal derivatives Z — N+. In this case Zq — N+. We integrate, again taking 
the lower limit to be the larger root of ^q — H"^ — 0, to get 



It will be shown below that this formula is also in resonable agreement with 
the analytical solution. 

The above approximations are subtle: we replace a cubic by a quadratic, chang- 
ing the integrand, but we also change the lower limit. These effects tend to 
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compensate, resulting in surprisingly accurate approximations. Moreover, it 
is found that the two approximations (71) and (72) have errors which are of 
opposite sign and approximately equal. Choosing Zq = {N + aJ)/2 in (67), 
we get the approximation 



By numerical experiment, we deduce the optimal value a — 0.458. Numerical 
results using the various approximations will be presented in the following 
section and (73) will be found to yield remarkably accurate results. 



6 Numerical Experiments 

We first compare the precession angle calculated using the exact analytical 
expression (57) with values extracted from a numerical integration of the three- 
wave equations (9)-(ll). For given N and J, the maximum value of the cubic 
$o(^) is at Zinax = |[2Ar - ^/WT^]. Thus, the maximum value of H is 



The three-wave equations were solved for a range of values < J < 1 and H 
covering the accessible parameter domain < H < Hq. Wc take = 1 in 
all cases; this is no loss of generality, as it is equivalent to a rescaling of the 
amplitudes by N'^^^ and of the time by N^^^. Prom the numerical solution, 
the major and minor axes 

^maj = 1^1 + 1^1 and A^in = \A\ - \B\ 

of the osculating or instantaneous ellipse (see [9]) were calculated as functions 
of time, and the precession angle was evaluated as the change in (p between 
successive maxima of Amaj- The precession angle was computed as a function 
of J and H. The results are presented in Fig. 5. The heavy line is Hq{J). 
The left-hand panel shows A0 calculated using the analytical formula (57). 
The precession angle vanishes for J = and is equal to 90° for H = 0. The 
center panel shows the angle calculated from the numerical solution of the 
three- wave equations. It is very similar to the analytical result. The right- 
hand panel shows the difference between the precession angle calculated from 
the numerical solution and the analytical formula. The values are generally 
very small (the contour interval in Fig. 5(C) is 0.1°). The maximum difference 
is 0.6° and the discrepancy may be ascribed to numerical noise. 




(73) 
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(A) A(|){ANA) (B) Ao{NUM) (C) |NUM-ANA] 




0.5 1 0.5 1 0.5 1 

J J J 



Fig. 5. Left-hand panel: precession angle calculated using the exact analyti- 
cal formula (57). Center panel: calculated from numerical integration of the 
three- wave equations. Note that = for J = and Acp = 90° for H = 0. 
Right-hand panel: difference in precession angle between the numerical and the 
analytical solution. 

6.1 Determination of the Precession Angle 



We now show that the envelope of the motion may be determined to high 
accuracy by using approximate formulas involving only elementary functions. 
We use the analytical values as a reference to evaluate the accuracy of the 
approximate formulas. In Fig. 6 the differences between the exact and ap- 
proximate expressions for A0 are shown. The absolute values of these errors 
are plotted. The maximum error in the DGC formula (Fig. 6(A)) is about 2.2°, 
and occurs for J fa | and H at its maximum permissible value. The error in 
the alternative formula (72) is of comparable magnitude, with a maximum of 
about 2.5° (Fig. 6(B)), but is of opposite sign. The optimal value a = 0.458 of 
the parameter in the formula (73) was found by experiment. Fig. 6(C) shows 
that this formula is significantly more accurate, with a maximum error less 
than 0.4°. This is a remarkable level of precision, considering the simplicity of 
the formula. The compensation of errors leads to what might be described as 
the unreasonable effectiveness of the approximation. 



6.2 Determination of the Pulsation Amplitude 



The extent to which energy is exchanged between the elastic and pendular 
modes of oscillation may be measured by the relative pulsation amplitude 
defined as 
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(A) |DGC-ANA| (B) |ALT-ANA| (C) |OPT - ANA| 
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Fig. 6. Differences in precession angle between tfiree approximate formulas and 

the analytical solution (57). (A) The DGC formula (71); (B) The alternative formula 
(72); (C) The optimum formula (73). Absolute values are shown. The signs of the 
errors of (71) and (72) are opposite. The contour interval is 0.1° in all panels. 

This quantity varies from P = for no energy exchange to P = 1 for maximal 
exchange. For H — 0, it reduces to P = 1 — J/N. Given the invariants A^, 
H and J, we may compute P by solving the cubic equation ^{Z) — where, 
as before, $(Z) = %{Z) - H^, with Z = |Cp and % defined by (26). For 
determination of the envelope, (74) is ideal. However, for the inverse problem, 
it must be simplified. Noting that Cf + C2 + = N, we may write the 
pulsation amplitude as 

N 



We have already introduced in (67) a quadratic ^0 which approximates the 
cubic $0 in the range [C|, Cj]. If we use the roots of ^0 — H'^ = as estimates 
of Cf and C|, an approximate expression for P may be obtained: 




(75) 



For fixed P this represents an eUipse in (J, i7)-space. The great advantage of 
(75) is that it can be used to solve the inverse problem. Two special cases 
follow immediately: when H = then P = 1 — J/iV (which is exact); when 
J = then P = 1 - 2H/N^/% (which is not exact). 

We plot the exact values of the pulsation amplitude, obtained by solving the 
cubic equation $(Z) = 0, in Fig. 7(A). Note that P = when H ^ Ho 
and P = 1 when H = J = 0. The approximate values calculated using 
(75) are shown in Fig. 7(B) and the difference (Papprox — -Pexact) in Fig- 7(C). 
The approximation is quite accurate when P is large. This is the region of 
primary physical interest, corresponding to strongly pulsating motion. For 
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(A) P-exact (B) P-approx (C) P-diff. 




0.5 1 0.5 1 0.5 1 

J J J 



Fig. 7. Pulsation amplitude. (A) P based on solving the cubic equation and using 
(74). (B) P from approximation (75). (C) Magnitude of difference between exact 
and approximate values The heavy curve is Hq{J). The contour interval is 0.1° in 
all panels. 

large values of H, the approximation is no longer valid. We have derived 
several other approximate expressions for P, which are more accurate, but 
also more complicated, than (75). 



6.3 Control of the envelope dynamics 



The approximate formulas allow us to control the pulsation and precession 
by a judicious choice of initial conditions. Recall that the precession angle is 
given, to high accuracy, by (73), which we write 

tanA0 = ^ (76) 



where Zq = [N + aJ)/2. This may be used in (75) to eliminate either H or 
J, yielding the two equations 

J IH 
P — 1 — — CSC and P — 1 ■= sec Ad . 



But these are instantly invertible, to give equations for J and H in terms of 
P and A0: 



J^N{1-P) sinA0 



and 



v^Ar(l -P)cosA(/). 



(77) 



To illustrate the effectiveness of these formulas, six values of the precession 
angle were chosen: Acj) e {10°, 20°, 30°, 45°, 60°, 90°}. We set iV = 1 and fixed 
the value of the pulsation amphtude to be P = 0.9. We then calculated J and 
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270 270 270 

Fig. 8. Polar plots of Amaj = |^| + \B\ against computed from the numerical 
solution of (9)-(ll) for six sets of initial conditions. For all cases, A'^ = 1 and 
P = 0.9, and J and H arc computed from (77). Top panels: G {10°, 20°, 30°}, 
bottom panels: Acj) G {45°, 60°, 90°}. The integration time in each case corresponds 
to a total precession of approximately 180°, and both Amaj and — Amaj are plotted. 

H from (77) and computed the numerical solution of the three-wave equations 
(9)-(ll). The initial value of \C\^ was taken to be the root C| of $(|C|2) = 
having intermediate algebraic magnitude. Then |Ap and \B\^ were obtained 
from the Manley-Rowe relations. The initial phases were all set to zero. Polar 
plots of ^maj against are shown in Fig. 8 (the integration time in each case 
corresponds to a total precession of about 180°, and both A^iaj and — Amaj are 
plotted). These plots represent the outer envelope of the horizontal projection 
of the trajectory of the pendulum bob. It is clear that the precession for the 
numerical solution is, in each case, close to the value used in (77). We also 
calculated the pulsation amplitude of the numerical solution and it was, in all 
cases, within 2% of the prescribed value P = 0.9. This confirms the effective- 
ness of the inversion formulas as a means of pre-determining the envelope of 
the motion. 

We note that, in general, the horizontal projection of the trajectory is not a 
closed curve, but densely covers a region of phase-space. The motion is not 
periodic but quasi-periodic. The horizontal projection is a closed curve only 
in the exceptional cases when A0 and 27r are commensurate, that is, when 
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QUALITATIVE ENVELOPE DYNAMICS 




Fig. 9. Qualitative features of the envelope dynamics of the swinging spring. 

their ratio is a rational number. In this case the motion is periodic and the 
horizontal footprint is a star-like graph, as illustrated in Fig. 8. The number 
of points in the star is the denominator of the rational number A0/27r {e.g., 
A4> = 80° yields a nine-pointed star). 



7 Conclusion 

We have presented a complete analytical solution of the three-wave equa- 
tions, which govern the small-amplitude dynamics of the resonant swinging 
spring. The periodic variation of the amplitudes is associated with the char- 
acteristic pulsation and precession of the system. Several analytical formulas 
for the precession angle have been presented. We have also derived simplified 
approximate expressions in terms of elementary functions. The optimal ap- 
proximation (73) has been shown by numerical experiments to be remarkably 
accurate, with a maximum error of only 0.4°. The amplitude of the pulsation 
envelope is determined from the roots of a cubic equation whose coefficients 
are defined by the invariants. Thus, we have provided a complete and positive 
answer to Question 1 posed in the Introduction. 

The inverse question. Question 2 in §1, has also been answered affirmatively. 
The approximate formulas (77) give values of J and H which lead to a solution 
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having the prescribed pulsation amphtude and precession angle. They are of 
high accuracy for strongly pulsating motion, which is the case of primary 
physical interest. 

The qualitative features of the envelope dynamics of the swinging spring are 
depicted schematically in Fig. 9. The axes are normalized angular momentum 
J/N and normalized Hamiltonian H/Hqq. The physically accessible domain 
is shaded. The bounding curve is H = Ho{J,N). The pulsation amplitude 
vanishes on this curve and the solutions are the elliptic-parabolic modes [13]. 
Regions of the parameter space are indicated where the pulsation amplitude 
and precession angle take large or small values. The corners of the accessible 
region represent special solutions. Thus, (J, H) = (0, -f^oo) corresponds to the 
cup-like and cap-like solutions of Vitt and Gorelik [17]. For (J, H) = {N, 0), the 
motion of the spring traces out a cone. Finally {J,H) = (0, 0) represents the 
homoclinic orbit, and includes the case of (unstable) pure vertical oscillations. 
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A The Nahm Equations and the Three-wave Equations 

The Nahm equations arc a set of integrablc equations for a three-vector of 
skew-Hermetian n x n matrices {Ti{s),T2{s),Ts{s)): 

^T, = [T„ n] = TjTk - nTj (A.l) 

where {ij k) is a cyclic permutation of (12 3). In the simplest case s e (—1, 1) 
and the matrices have simple poles as s — > ±1. 

The Nahm equations were originally discovered because it is possible to use 
solutions to the equations to construct solutions of the Bogomolny equation 

[16]. These solutions are called Bogomolny-Prasad-Sommcrficld monopolcs. 
The Bogomolny equation occurs as a super-symmetry or minimum energy 
condition in Yang- Mills Higgs theory and is of interest to theoretical particle 
physicists. 



28 



There is a Lax formulation of the Nahm equations and an associated Lax 
curve of genus (n — 1)^. The n = 2 case is elhptic and the sohitions arc clhptic 
functions; in fact, for n = 2 the Nahm equations reduce to the Euler-Poinsot 
equations and are easily solved. Surprisingly, it is sometimes also the case that 
the Nahm equations forn > 2 can be solved in terms of elliptic functions. This 
happens when the solution has a symmetry and the quotient of the Lax curve 
by that symmetry gives a genus-one surface. These symmetries of the Nahm 
matrices correspond to spatial symmetries of the corresponding monopoles. 
The group elements act both by conjugation on the Nahm matrices and by 
rotation of the three- vector of matrices [8,11]. 



One example is n 
matrices to 



3 D2 symmetry [10]. The symmetry reduces the Nahm 



^ F* ^ 

Fi Fi 
yo F* ) 



1 

V2 



-F* F* 
-Fo 



^ -i^iFs} -^{F:,}^ 
T3= 

^ Q{Fs} iU{Fs} J 

Substituting these matrices into the Nahm equations gives 
dFi 



ds 



f;f* 



(A.2) 



and two others by cyclic permutation. These equations are the 'explosive inter- 
action' three- wave equations identified in [2] . They are related to the equations 
studied in the present paper hy s — it with Fi — A*, F2 — B* and F3 = C. 



B Relationship between Jacobi and Weierstrass forms of the Pre- 
cession Angle. 



To relate the expression (57) obtained by means of Jacobi's elhptic functions 
to the expression (48) in terms of the Weierstrass form, we introduce auxiliary 
constants 0?+ and ci_ defined by 
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Note that since 7I. > k'^, these constants are complex (d+ and d- lie on the 
hne between K and K + iK'). It follows that 



with similar expressions involving d-. The first of (55) may be written 



It may be shown without difficulty, using equation (50), that 
= -\-ik^ sn d+ cn 0?+ dn 0?+ . 

Then the solution (56) for ^ may be written 

e-eo = -A+s-mi(s,d+,A;). (B.l) 

where ni(s, (i+,/c) is another standard form (Jacobi's form) for the elliptic 
integral of the third kind ([19], §22.74): 

k"^ sn (ii cn dn d< sn^s 



, . /■ fc sn d+ cn d+ dn ct+ sn s 

^) = i 1-fe^sn^d.sn^. ■ 



The eUiptic integral of the second kind is defined (with x — sn.z) as 



dv?zdz = j sl-^—^dx. 







The complete integral is denoted E = E{K). E{z) is not periodic; the periodic 
component is represented by Jacobi's zeta function 

Z{z) = E{z)-Ez/K . 



This is an odd function with period 2K. It is related to Jacobi's theta function, 
also having period 2K, by 

z{z) ^ j-^\oge{z) . 
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The elliptic integral of the third kind may now be expressed as follows: 

ni(^, a,k)^l log + zZ{a) (B.2) 



For the complete form of the integral, when z = K, the logarithmic term 
vanishes: 

U,(a,k)^KZ(a)^K^, (B.3) 



Using this in (B.l), we obtain the change over a half-period K: 
|AC = -KX+ - iKZ{d+) . 

Finally, using the analogous expression for Arj, we get the precession angle 

(B.4) 



A(j) = -K{X+ - A_) - iK{Z{d+) - Z{dJ)) 



This is the change over the period 2K (for s) or 2K/v^/ N (for t). The struc- 
tural similarity between this expression and the result (48) in terms of Weier- 
strass functions is immediate. 

When J = 0, we have A_)_ = A_ and 7+ = 7_, so that 0?+ = d_ and (B.4) 
implies A0 = 0. For H — we have \+ — — 0, d+ — K and d_ = K + iK'. 
Then using the relation 

, ivr cnMdn?i 
with u = K, it follows immediately that A0 = 7r/2, in agreement with (49). 



C Approximation Integral with Best Fit at Z = N_|_. 

We approximate $0 — Z{Z — N^){Z — N^) by a quadratic with a root at 
Z-N+: 

^ Zo{Z - N+){Z - Zi) (C.l) 
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To obtain the best fit at Z = N^, we choose Zq and Zi so that 

^o{Z) - MZ) = 0({Z- iV+)3) . (C.2) 



This imphes Zq = + J and Zi = A^^/(A^+ + J). Using the software package 
Maple, it is possible to evaluate the resulting approximate integral 

A^^j^/ „... (C.3, 



2 



where C+ is the larger root of \E'o — = 0. The result is 



4y/j^~^lP V 2V2VJ^ - HWN + 3J 

This gives a rather good approximation: the maximum error is about 1.2°. 
However, it is not a convenient expression because the factor \/ — H"^ is 
sometimes real and sometimes imaginary. Moreover, the expression cannot 
easily be inverted to give if or J in terms of A0. In fact, unless Zi = A^_, 
any approximating quadratic (C.l) will give an expression with this problem. 
Thus, in §5, we choose Zi = N_. 
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